close all
clear
clc
%%%%%%%%%%%%%%%%%%%%%
%     Parameter     %
%%%%%%%%%%%%%%%%%%%%%
a=0.2;                      %lowpass filter not used
slideWindow=15;             %delta function forwarding and tailing window
cutoff=10;                  %peak selection cutoff critias
excludePercentage=0.1;      %fitting excluding percentage
fitWindow=15;
cutoffForce=8.01;
stepCutoffMagnet=8.01;
internalTTest=1;            %Single T test
internalTTestAlpha=0.01;
position_overstretching=12.41;
aa_length=0.5;    % the pesistent length of a amino acid
filelist=importdata('G:\unfolding_vinculin\filelistnew1.txt');
only_unfolding=1;% only process the data in the unfolding part of spectrum
with_smoothing=0;%whether to smooth rawdata
drift_correction=0; %whether to correct linear drift in long time measurements

%note the absolute extension of individual trajectories at specified
%magnetic positions
output_mean_extension=1;


for currentfile = filelist'
    
    name=char(currentfile);
    
    
    %% Data processing
    steps_data=[];
    filename=sprintf('%s%s',name,'.txt')
    output_filename=sprintf('%s%s',name,'-out.txt');
    
    outfid=fopen(output_filename,'w');
    %c=65/(exp(-(13.5-position_overstretching)/0.36)+0.48*exp(-(13.5-position_overstretching)/1.12));
    c=280;%for M270
    datastruct=importdata(filename);
    dataraw=datastruct.data;
    egxT=dataraw(:,1)';
    egMagneticT=dataraw(:,2);
    egRawDynT=dataraw(:,3);
    %plot(egxT); % Check equal spacing
    
    %Mingxi data: A. linear drift remover
    %             B. force altering regions
    %             C. points for fitting selection
    
    if drift_correction
        linearFitP=polyfit(egxT(egMagneticT==min(egMagneticT))',egRawDynT(egMagneticT==min(egMagneticT)),1);
        egfittedDriftT=polyval(linearFitP, egxT)';
        disp('Determining the drift\n');
        linearFitP % display the drift
        egRawDynTShift=egRawDynT-egfittedDriftT;
    else
        egRawDynTShift=egRawDynT;
    end
    %figure;
    %plot(egxT,egRawDynT);
    %hold on; plot(egxT,egfittedDriftT,'r','LineWidth',3);
    %figure;
    %plot(egxT,egRawDynT-egfittedDriftT,'b');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %cutoffForce=min(egMagneticT);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %process if only unfolding signal is processed
    tsig=egMagneticT>cutoffForce;
    dsig=diff([0;tsig;0]);
    startIndex=find(dsig > 0);
    endIndex=find(dsig < 0)-2;
    
    
    hFig = figure();
    set(hFig, 'Visible', 'on')
    
    
    for forceAlteringRegionIndex=1:length(startIndex)
        
       
        region_end_index=endIndex(forceAlteringRegionIndex);
        region_start_index=startIndex(forceAlteringRegionIndex);
        if only_unfolding    %with this option only the signals in the unfolding region of the trajectory is analyzed
            region_Magnetic=egMagneticT(startIndex(forceAlteringRegionIndex):endIndex(forceAlteringRegionIndex)); % find the magnetic positions in the region
            region_magsig=diff([0;region_Magnetic;0]);
            if ~isempty(find(region_magsig<0,1,'first'))
                region_end_index=region_start_index+find(region_magsig<0,1,'first')-1;
            end
        end
        %region_start_index
        %region_end_index
        egx=egxT(region_start_index:region_end_index);
        egRawDyn=egRawDynTShift(region_start_index:region_end_index);
        %scaled magnetic Location
        egMagnetic=egMagneticT(region_start_index:region_end_index);
        egForce=c*(exp(-(13.5-egMagnetic)./0.36)+0.48.*exp(-(13.5-egMagnetic)./1.12));
        
        %egTruDyn=(egMagnetic-min(egMagnetic))*(max(egRawDyn)-min(egRawDyn))/(max(egMagnetic)-min(egMagnetic));
        %figure;
        %plot(egxT,egRawDynTShift);
        %hold on;plot(egx,egRawDyn,'r');
        
        
        %% Assign
        x=egx;
        RawDyn=egRawDyn;
        %TruDyn=egTruDyn;
        %% Filtering By lowpass filter (Freqency**)
        % Without filter seems better for loading rate data (or change parameter)
        FilDynA=filter(a,[1 a-1],RawDyn);
        FilDynB=smooth(RawDyn,slideWindow,'sgolay'); % SG smooth
        %Average the extension for individual forces in the unfolding
        %trajectory
        [Fa,Fb,Fc] =unique(egForce);
        extension_forcemean=[Fa,accumarray(Fc,egRawDyn,[],@mean)];
        figure(1)
        plot(x,RawDyn,'b',x,FilDynA,'r');
        hold on; plot(x,FilDynB,'g');
        figure(2)
        plot(extension_forcemean(:,1),extension_forcemean(:,2),'b-', 'Color', rand(1,3));
        hold on;
       %plot the average force extension curve
        
       %keep track of the extension 
       outputmatrix=[repmat(forceAlteringRegionIndex,length(extension_forcemean),1) extension_forcemean];
        if output_mean_extension
            fprintf(outfid,'%d\t%f\t%f\n',outputmatrix');
        end

        
 
    end
    figure_filename = sprintf('%s%s',name,'.png');
   
%    print(gcf, '-dpng','-opengl', figure_filename, '-r640');
    
   fclose(outfid);
    
    
end
